# Carga de librerías
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.2 ✓ purrr 0.3.4
## ✓ tibble 3.0.3 ✓ dplyr 1.0.0
## ✓ tidyr 1.1.0 ✓ stringr 1.4.0
## ✓ readr 1.3.1 ✓ forcats 0.5.0
## ── Conflicts ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
library(tidymodels)
## ── Attaching packages ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidymodels 0.1.1 ──
## ✓ broom 0.7.0 ✓ recipes 0.1.13
## ✓ dials 0.0.9 ✓ rsample 0.0.8
## ✓ infer 0.5.3 ✓ tune 0.1.1
## ✓ modeldata 0.1.0 ✓ workflows 0.2.1
## ✓ parsnip 0.1.4 ✓ yardstick 0.0.7
## ── Conflicts ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidymodels_conflicts() ──
## x scales::discard() masks purrr::discard()
## x dplyr::filter() masks stats::filter()
## x recipes::fixed() masks stringr::fixed()
## x dplyr::lag() masks stats::lag()
## x yardstick::spec() masks readr::spec()
## x recipes::step() masks stats::step()
library(GGally)
## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2
library(corrr)
library(knitr)
library(ggrepel)
library(kableExtra)
##
## Attaching package: 'kableExtra'
## The following object is masked from 'package:dplyr':
##
## group_rows
library(here)
## here() starts at /home/daro/cdatos/TPs/EEA/TP2
El dataset de precios de inmuebles proviene de Properati. El mismo ya fue filtrado por los docentes y a su vez se encuentra particionado en subconjuntos de training y test.
train_ds <- read_csv(here("/ds/ar_properties_train.csv"))
## Parsed with column specification:
## cols(
## id = col_character(),
## l3 = col_character(),
## rooms = col_double(),
## bathrooms = col_double(),
## surface_total = col_double(),
## surface_covered = col_double(),
## price = col_double(),
## property_type = col_character()
## )
head(train_ds)
## # A tibble: 6 x 8
## id l3 rooms bathrooms surface_total surface_covered price property_type
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 r22O… Alma… 1 1 40 37 92294 Departamento
## 2 atZQ… Alma… 1 1 49 44 115000 Departamento
## 3 Tur/… Alma… 1 1 40 37 88900 Departamento
## 4 PNtN… Alma… 1 1 40 37 88798 Departamento
## 5 Lhg4… Alma… 1 1 49 44 110975 Departamento
## 6 oeLh… Alma… 1 1 40 37 99000 Departamento
dim_desc(train_ds)
## [1] "[32,132 x 8]"
summary(train_ds)
## id l3 rooms bathrooms
## Length:32132 Length:32132 Min. :1.000 Min. :1.000
## Class :character Class :character 1st Qu.:2.000 1st Qu.:1.000
## Mode :character Mode :character Median :3.000 Median :1.000
## Mean :2.722 Mean :1.427
## 3rd Qu.:3.000 3rd Qu.:2.000
## Max. :8.000 Max. :5.000
## surface_total surface_covered price property_type
## Min. : 28 Min. : 26.00 Min. : 69500 Length:32132
## 1st Qu.: 46 1st Qu.: 42.00 1st Qu.:120000 Class :character
## Median : 65 Median : 58.00 Median :169000 Mode :character
## Mean : 80 Mean : 70.61 Mean :214991
## 3rd Qu.: 99 3rd Qu.: 85.00 3rd Qu.:260000
## Max. :320 Max. :265.00 Max. :950000
unique(train_ds$l3)
## [1] "Almagro" "Palermo" "Caballito"
## [4] "Villa Crespo" "Villa Urquiza" "Barracas"
## [7] "Nuñez" "Belgrano" "Floresta"
## [10] "Recoleta" "Colegiales" "Parque Chas"
## [13] "Barrio Norte" "Villa Devoto" "Villa Ortuzar"
## [16] "Velez Sarsfield" "Parque Chacabuco" "Villa Pueyrredón"
## [19] "Villa Real" "Once" "Flores"
## [22] "San Telmo" "Las Cañitas" "Villa Santa Rita"
## [25] "Villa del Parque" "Retiro" "Parque Centenario"
## [28] "Chacarita" "Abasto" "San Cristobal"
## [31] "Liniers" "Balvanera" "Boedo"
## [34] "Villa General Mitre" "Saavedra" "Parque Patricios"
## [37] "Boca" "Paternal" "Monserrat"
## [40] "San Nicolás" "Parque Avellaneda" "Centro / Microcentro"
## [43] "Puerto Madero" "Congreso" "Monte Castro"
## [46] "Mataderos" "Coghlan" "Versalles"
## [49] "Villa Lugano" "Catalinas" "Villa Luro"
## [52] "Constitución" "Pompeya" "Tribunales"
## [55] "Agronomía" "Villa Soldati" "Villa Riachuelo"
unique(train_ds$property_type)
## [1] "Departamento" "PH" "Casa"
apply(train_ds, 2, function(x) any(is.na(x)))
## id l3 rooms bathrooms surface_total
## FALSE FALSE FALSE FALSE FALSE
## surface_covered price property_type
## FALSE FALSE FALSE
table(train_ds$property_type)
##
## Casa Departamento PH
## 809 28203 3120
tail(table(train_ds$l3)[], 10)
##
## Villa General Mitre Villa Lugano Villa Luro Villa Ortuzar
## 99 113 264 109
## Villa Pueyrredón Villa Real Villa Riachuelo Villa Santa Rita
## 285 55 13 128
## Villa Soldati Villa Urquiza
## 9 1319
Vamos a generar un modelo que utilice todas las covariables presentes en nuestro dataset de training para buscar predecir el precio.
modelo_todas_covariables <- lm(
price ~ rooms + bathrooms + surface_total + surface_covered + property_type + l3,
data = train_ds
)
tidy_mtc <- tidy(modelo_todas_covariables, conf.int = TRUE)
tidy_mtc
## # A tibble: 63 x 7
## term estimate std.error statistic p.value conf.low conf.high
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) -109144. 5784. -18.9 5.28e- 79 -120480. -97808.
## 2 rooms -4360. 530. -8.23 2.01e- 16 -5399. -3321.
## 3 bathrooms 34367. 768. 44.7 0. 32862. 35873.
## 4 surface_total 891. 28.2 31.5 3.84e-215 835. 946.
## 5 surface_covered 1498. 34.3 43.6 0. 1431. 1565.
## 6 property_typeDepar… 91486. 2660. 34.4 1.35e-254 86272. 96700.
## 7 property_typePH 46220. 2752. 16.8 5.06e- 63 40826. 51615.
## 8 l3Agronomía 1241. 10652. 0.117 9.07e- 1 -19638. 22120.
## 9 l3Almagro -3480. 5182. -0.672 5.02e- 1 -13636. 6676.
## 10 l3Balvanera -25835. 5491. -4.70 2.55e- 6 -36598. -15072.
## # … with 53 more rows
El valor de la ordenada al origen (-109.143,78) es el valor del precio esperado para una propiedad sin ambientes o superficie y que no pertenece a ningún tipo (Casa, PH o Departamento) y que no se encuentra en ninguna zona (es una inferencia más teórica de nuestro modelo dado que no es representativo en la realidad).
El coeficiente estimado de ambientes es -4359,73. Esto implica que manteniendo constantes todas nuestras otras covariables, un aumento de un ambiente a nuestra propiedad corresponde a una quita de $4359, en promedio en el precio de una propiedad.
El coeficiente estimado de baños es 34.367,21. Esto implica que manteniendo constantes todas nuestras otras covariables, un aumento de un baño a nuestra propiedad corresponde a una aumento de $34.367, en promedio en el precio de una propiedad.
El coeficiente estimado de superficie total es 890,56. Esto implica que manteniendo constantes todas nuestras otras covariables, un aumento en una unidad de superficie (m2) a nuestra propiedad corresponde a una aumento de $890, en promedio en el precio de una propiedad.
El coeficiente estimado de superficie cubierta es 1.497,93. Esto implica que manteniendo constantes todas nuestras otras covariables, un aumento en una unidad de superficie cubierta (m2) a nuestra propiedad corresponde a una aumento de $1.498, en promedio en el precio de una propiedad.
Los coeficientes de tipo de propiedad (Departamento y PH) indican cuanto aumenta la función de respuesta (precio) para un departamento o un ph en comparación a una casa (categoría basal). Sus valores de β respectivos son para los departamentos 91.485,94 y para los PHs 46.220,04
Finalmente tenemos 56 coeficientes para la variable l3 que representan en cuanto aumenta o se reduce el precio de una propiedad dependiendo de la zona específica a la cual pertenece.
Observamos en nuestro modelo dos tipos distintos de variables dummies. Aquellas que responden a la variable property_type y aquellas que responden a la variable l3.
Podemos observar que los p-valores asociados a los coeficientes de la variable property_type, resultan estadísticamente significativos. En cuanto a aquellos p-valores asociados a los coeficientes de la variable l3 nos encontramos con que en algunos casos resultan estadísticamente significativos (l3Balvanera con un p-valor de 0,000002) pero en otros casos no lo es (l3Agronomía con un p-valor de 0,9).
Además podemos observar que los intervalos de confianza para las variables dummy que responden a property_type no contienen al 0 y en lineas generales, aquellos que responden a la variable l3 si lo contienen.
Vamos a aplicar un Test F (test de ignificatividad global) para medir la significatividad conjunta de las variables categóricas en nuestro modelo.
tidy(
anova(modelo_todas_covariables)
)
## # A tibble: 7 x 6
## term df sumsq meansq statistic p.value
## <chr> <int> <dbl> <dbl> <dbl> <dbl>
## 1 rooms 1 2.18e14 2.18e14 48554. 0
## 2 bathrooms 1 1.09e14 1.09e14 24391. 0
## 3 surface_total 1 7.74e13 7.74e13 17276. 0
## 4 surface_covered 1 1.73e13 1.73e13 3863. 0
## 5 property_type 2 2.03e13 1.01e13 2265. 0
## 6 l3 56 5.90e13 1.05e12 235. 0
## 7 Residuals 32069 1.44e14 4.48e 9 NA NA
La tabla de ANOVA muestra que tanto la variable categorica property_type como l3 en su conjunto resultan estadísticamente significativas (p-valor < 0.05). Por ende buscaremos determinar cuales de las comparaciones entre grupos son estadísticamente significativas.
Graficamos los coeficientes estimados.
ggplot(
tidy_mtc,
aes(
estimate,
term,
xmin = conf.low,
xmax = conf.high,
height = 0
)
) +
geom_point(
color = "forestgreen",
size=2
) +
geom_vline(
xintercept = 0,
lty = 4,
color = "black"
) +
geom_errorbarh(
color = "forestgreen",
size=1
) +
theme_bw() +
labs(
y = "Coeficientes β",
x = "Estimación"
)
De analizar el grafico de los coeficientes podemos visualizar los intervalos de confianza de las distintas variables. Claramente podemos observar como hay muchas variables dummy asociadas a la variable l3 que incluyen el 0 dentro de sus intervalos de confianza y por ende podemos asumir que estas variables no son significativas para explicar el precio de las propiedades.
Y ahora analizamos las medidas del resumen global del modelo.
glance(modelo_todas_covariables)
## # A tibble: 1 x 12
## r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0.777 0.777 66933. 1803. 0 62 -4.03e5 8.05e5 8.06e5
## # … with 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>
Analizando el valor del \(Rˆ2\) ajustado del modelo, podemos ver que el mismo parece explicar un buen porcentaje de variabilidad del conjunto.
### Generación de Modelo sin l3
Ahora vamos a generar un modelo sin la variable l3 para comparar y analizar el impacto de esta variable en nuestro modelo previo.
modelo_sin_l3 <- lm(
price ~ rooms + bathrooms + surface_total + surface_covered + property_type,
data = train_ds
)
tidy_msl3 <- tidy(modelo_sin_l3, conf.int = TRUE)
tidy_msl3
## # A tibble: 7 x 7
## term estimate std.error statistic p.value conf.low conf.high
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) -128676. 3333. -38.6 1.07e-318 -135210. -122143.
## 2 rooms -13657. 617. -22.1 8.68e-108 -14867. -12448.
## 3 bathrooms 42564. 899. 47.3 0. 40801. 44327.
## 4 surface_total 856. 33.0 25.9 1.40e-146 791. 920.
## 5 surface_covered 1819. 40.0 45.5 0. 1740. 1897.
## 6 property_typeDepart… 133033. 3049. 43.6 0. 127056. 139010.
## 7 property_typePH 66457. 3233. 20.6 2.59e- 93 60121. 72793.
Podemos observar como en nuestro nuevo modelo, los p-valores de todos los coeficientes asociados a nuestras variables resultan estadísticamente significativos (p-valor < 0.05). Además ninguno de los intervalos de confianza de nuestros coeficientes incluyen el 0.
Vamos a aplicar un Test F (test de ignificatividad global) para medir la significatividad conjunta de las variables categóricas en nuestro modelo.
tidy(
anova(modelo_sin_l3)
)
## # A tibble: 6 x 6
## term df sumsq meansq statistic p.value
## <chr> <int> <dbl> <dbl> <dbl> <dbl>
## 1 rooms 1 2.18e14 2.18e14 34486. 0
## 2 bathrooms 1 1.09e14 1.09e14 17324. 0
## 3 surface_total 1 7.74e13 7.74e13 12271. 0
## 4 surface_covered 1 1.73e13 1.73e13 2744. 0
## 5 property_type 2 2.03e13 1.01e13 1609. 0
## 6 Residuals 32125 2.03e14 6.31e 9 NA NA
Los resultados del ANOVA muestran que todas las variables utilizadas son estadisticamente significativas (respaldando lo visto en el anterior modelo).
Y ahora analizamos comparativamente las medidas del resumen global de nuestros 2 modelos.
modelos = list(
modelo_todas_covariables = modelo_todas_covariables,
modelo_sin_l3 = modelo_sin_l3
)
purrr::map_df(modelos, broom::glance, .id = "model")
## # A tibble: 2 x 13
## model r.squared adj.r.squared sigma statistic p.value df logLik AIC
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 mode… 0.777 0.777 66933. 1803. 0 62 -4.03e5 8.05e5
## 2 mode… 0.686 0.686 79419. 11674. 0 6 -4.08e5 8.16e5
## # … with 4 more variables: BIC <dbl>, deviance <dbl>, df.residual <int>,
## # nobs <int>
Como se puede observar, el modelo_todas_covariables es el modelo que mejor explica la variabilidad en nuestro conjunto (evaluamos el coeficiente \(Rˆ2\) ajustado para ambos).
En nuestro analisis previo, observamos como algunos barrios (variable l3) eran significativos y otros no. Por ende, vamos a proponer la generación de una nueva variable barrios que nos permite agrupar mas eficientemente nuestros datos.
Para hacerla, vamos a generar otra variable precio_m2 que contiene el valor de price sobre el total de metros cuadrados de superficie (surface_total). Luego, sobre esta variable, agruparemos los datos en tres categorías (variable barrios): - precio_bajo: aquellas propiedades con valores iguales o inferiores al Q1 de nuestra variable precio_m2 - precio_medio: aquellas propiedades con valores superiores al Q1 y con valores inferiores o iguales al Q3 de nuestra variable precio_m2 - precio_alto: aquellas propiedades con valores superiores al Q3 de nuestra variable precio_m2
train_ds = train_ds %>%
mutate(
precio_m2 = price / surface_total,
barrios = factor(
case_when(
precio_m2 <= quantile(precio_m2)[2] ~ "precio_bajo",
precio_m2 > quantile(precio_m2)[2] & precio_m2 <= quantile(precio_m2)[4] ~ "precio_medio",
precio_m2 > quantile(precio_m2)[4] ~ "precio_alto"
)
)
)
train_ds
## # A tibble: 32,132 x 10
## id l3 rooms bathrooms surface_total surface_covered price
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 r22O… Alma… 1 1 40 37 92294
## 2 atZQ… Alma… 1 1 49 44 115000
## 3 Tur/… Alma… 1 1 40 37 88900
## 4 PNtN… Alma… 1 1 40 37 88798
## 5 Lhg4… Alma… 1 1 49 44 110975
## 6 oeLh… Alma… 1 1 40 37 99000
## 7 Fmtr… Alma… 1 1 40 37 96984
## 8 QohS… Pale… 1 1 40 34 99500
## 9 WEJo… Pale… 1 1 36 33 121600
## 10 ByY9… Pale… 3 2 90 90 285000
## # … with 32,122 more rows, and 3 more variables: property_type <chr>,
## # precio_m2 <dbl>, barrios <fct>
Ahora vamos a generar un nuevo modelo incluyendo nuestras nueva covariables barrios.
modelo_barrios <- lm(
price ~ rooms + bathrooms + surface_total + surface_covered + property_type + barrios ,
data = train_ds
)
tidy_mb <- tidy(modelo_barrios, conf.int = TRUE)
tidy_mb
## # A tibble: 9 x 7
## term estimate std.error statistic p.value conf.low conf.high
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) -3082. 2574. -1.20 2.31e- 1 -8127. 1962.
## 2 rooms -3831. 458. -8.37 6.08e- 17 -4728. -2934.
## 3 bathrooms 21781. 674. 32.3 2.58e-225 20459. 23102.
## 4 surface_total 1561. 24.9 62.8 0. 1512. 1610.
## 5 surface_covered 1103. 29.8 37.0 1.22e-293 1044. 1161.
## 6 property_typeDepart… 87739. 2266. 38.7 1.02e-320 83299. 92180.
## 7 property_typePH 59258. 2379. 24.9 1.14e-135 54596. 63921.
## 8 barriosprecio_bajo -165257. 1017. -163. 0. -167250. -163264.
## 9 barriosprecio_medio -93607. 810. -116. 0. -95195. -92020.
Podemos apreciar como los p-valores de los coeficientes resultan significativos. Sin embargo el valor de mi coeficiente \(B_0\) parece no ser significativo (p-valor > 0.05). Más aún, podemos comprobar que su intervalo de confianza incluye el 0.
Comparamos nuestros 3 modelos.
modelos = list(
modelo_todas_covariables = modelo_todas_covariables,
modelo_sin_l3 = modelo_sin_l3,
modelo_barrios = modelo_barrios
)
purrr::map_df(modelos, broom::glance, .id = "model") %>%
arrange(adj.r.squared)
## # A tibble: 3 x 13
## model r.squared adj.r.squared sigma statistic p.value df logLik AIC
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 mode… 0.686 0.686 79419. 11674. 0 6 -4.08e5 8.16e5
## 2 mode… 0.777 0.777 66933. 1803. 0 62 -4.03e5 8.05e5
## 3 mode… 0.830 0.830 58434. 19576. 0 8 -3.98e5 7.97e5
## # … with 4 more variables: BIC <dbl>, deviance <dbl>, df.residual <int>,
## # nobs <int>
Podemos observar como el modelo que mejor explica la variabilidad, \(Rˆ2\) ajustado, es nuestro modelo_barrios seguido en orden por modelo_todas_covariables y modelo_sin_l3.
Más aún, como el modelo_barrios tiene coeficientes que son estadísticamente significativos, resulta más útil para el análisis del precio de nuestras propiedades (a pesar que la información de los valores de l3 podría ser más util en un análisis mas profundo).
Dada la fuerte correlación entre las variables surface_total y surface_covered, vamos a generar una nueva variable llamada sup_descubierta que representa la diferencia entre la superficie total y la cubierta.
train_ds = train_ds %>%
mutate(sup_descubierta = surface_total - surface_covered)
train_ds
## # A tibble: 32,132 x 11
## id l3 rooms bathrooms surface_total surface_covered price
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 r22O… Alma… 1 1 40 37 92294
## 2 atZQ… Alma… 1 1 49 44 115000
## 3 Tur/… Alma… 1 1 40 37 88900
## 4 PNtN… Alma… 1 1 40 37 88798
## 5 Lhg4… Alma… 1 1 49 44 110975
## 6 oeLh… Alma… 1 1 40 37 99000
## 7 Fmtr… Alma… 1 1 40 37 96984
## 8 QohS… Pale… 1 1 40 34 99500
## 9 WEJo… Pale… 1 1 36 33 121600
## 10 ByY9… Pale… 3 2 90 90 285000
## # … with 32,122 more rows, and 4 more variables: property_type <chr>,
## # precio_m2 <dbl>, barrios <fct>, sup_descubierta <dbl>
Ahora con nuestra nueva variable, vamos a generar un nuevo modelo que incluya la variable sup_descubierta en reemplazo de la variable surface_total.
modelo_sup_descubierta <- lm(
price ~ rooms + bathrooms + sup_descubierta + surface_covered + property_type + precio_m2 + barrios ,
data = train_ds
)
tidy_msd <- tidy(modelo_sup_descubierta, conf.int = TRUE)
tidy_msd
## # A tibble: 10 x 7
## term estimate std.error statistic p.value conf.low conf.high
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) -343392. 2485. -138. 0. -3.48e5 -338521.
## 2 rooms 1064. 313. 3.41 6.61e- 4 4.52e2 1677.
## 3 bathrooms 12549. 461. 27.2 3.96e-161 1.16e4 13453.
## 4 sup_descubierta 1970. 17.1 115. 0. 1.94e3 2004.
## 5 surface_covered 2611. 10.5 249. 0. 2.59e3 2632.
## 6 property_typeDepa… 68762. 1545. 44.5 0. 6.57e4 71790.
## 7 property_typePH 51785. 1619. 32.0 6.01e-221 4.86e4 54959.
## 8 precio_m2 92.2 0.478 193. 0. 9.13e1 93.2
## 9 barriosprecio_bajo 19276. 1180. 16.3 1.00e- 59 1.70e4 21590.
## 10 barriosprecio_med… 18396. 800. 23.0 5.87e-116 1.68e4 19965.
Podemos observar como los p-valores de todos nuestros coeficientes resultan significativos en nuestro nuevo modelo.
tidy(
anova(modelo_sup_descubierta)
)
## # A tibble: 8 x 6
## term df sumsq meansq statistic p.value
## <chr> <int> <dbl> <dbl> <dbl> <dbl>
## 1 rooms 1 2.18e14 2.18e14 137576. 0.
## 2 bathrooms 1 1.09e14 1.09e14 69112. 0.
## 3 sup_descubierta 1 3.81e12 3.81e12 2407. 0.
## 4 surface_covered 1 9.09e13 9.09e13 57490. 0.
## 5 property_type 2 2.03e13 1.01e13 6417. 0.
## 6 precio_m2 1 1.51e14 1.51e14 95488. 0.
## 7 barrios 2 8.63e11 4.32e11 273. 2.78e-118
## 8 Residuals 32122 5.08e13 1.58e 9 NA NA
Y al ejecutar el Test F, validamos el mismo.
Ahora comparamos todos los modelos que generamos.
modelos = list(
modelo_todas_covariables = modelo_todas_covariables,
modelo_sin_l3 = modelo_sin_l3,
modelo_barrios = modelo_barrios,
modelo_sup_descubierta = modelo_sup_descubierta
)
purrr::map_df(modelos, broom::glance, .id = "model") %>%
arrange(adj.r.squared)
## # A tibble: 4 x 13
## model r.squared adj.r.squared sigma statistic p.value df logLik AIC
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 mode… 0.686 0.686 79419. 11674. 0 6 -4.08e5 8.16e5
## 2 mode… 0.777 0.777 66933. 1803. 0 62 -4.03e5 8.05e5
## 3 mode… 0.830 0.830 58434. 19576. 0 8 -3.98e5 7.97e5
## 4 mode… 0.921 0.921 39763. 41717. 0 9 -3.86e5 7.72e5
## # … with 4 more variables: BIC <dbl>, deviance <dbl>, df.residual <int>,
## # nobs <int>
Podemos observar como nuestro nuevo modelo modelo_sup_descubierta es el que mejor explica la variabilidad.
Vamos a evaluar nuestro modelo lineal con la variable sup_descubierta
plot(modelo_sup_descubierta)
- Residuos vs valores predichos: parece existir una cierta estructura en los datos. La linealidad parece que no se respeta. Vemos la existencía de heterocedasticidad a medida que aumentan las predicciones. Se resaltan 2 observaciones. - Normal QQ plot: los extremos no se ajustan a la distribución teórica, tanto inferior izquierdo como superior derecho. - Residual vs leverage:Existen cuatro puntos con un leverage bastante alto. Y se identifican 2 puntos como posibles outliers. .
Vamos a generar un modelo para predecir el logaritmo natural de la variable price \[ log(price) = \beta_0 + \beta_1log(rooms) + \beta_2log(bathrooms) + \beta_3log(surface\_covered) + \beta_4property\_type + \beta_5barrio + \beta_6surface\_patio \]
Para esto vamos a generar nuestros nuevas variables
modelo_log = lm(
log(price) ~ log(rooms) + log(bathrooms) + log(surface_covered) + property_type + barrios + sup_descubierta,
data = train_ds
)
tidy_ml <- tidy(modelo_log, conf.int = TRUE)
tidy_ml
## # A tibble: 9 x 7
## term estimate std.error statistic p.value conf.low conf.high
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 8.53 0.0155 549. 0. 8.50 8.56
## 2 log(rooms) -0.0211 0.00315 -6.68 2.36e- 11 -0.0273 -0.0149
## 3 log(bathrooms) 0.0694 0.00324 21.4 2.68e-101 0.0630 0.0757
## 4 log(surface_covered) 0.901 0.00371 243. 0. 0.894 0.908
## 5 property_typeDepart… 0.171 0.00611 27.9 1.00e-169 0.159 0.183
## 6 property_typePH 0.0921 0.00644 14.3 3.61e- 46 0.0794 0.105
## 7 barriosprecio_bajo -0.708 0.00278 -254. 0. -0.713 -0.702
## 8 barriosprecio_medio -0.359 0.00222 -162. 0. -0.364 -0.355
## 9 sup_descubierta 0.00694 0.0000682 102. 0. 0.00680 0.00707
tidy(
anova(modelo_log)
)
## # A tibble: 7 x 6
## term df sumsq meansq statistic p.value
## <chr> <int> <dbl> <dbl> <dbl> <dbl>
## 1 log(rooms) 1 3724. 3724. 145610. 0
## 2 log(bathrooms) 1 1550. 1550. 60606. 0
## 3 log(surface_covered) 1 1381. 1381. 54013. 0
## 4 property_type 2 172. 86.2 3371. 0
## 5 barrios 2 1481. 740. 28948. 0
## 6 sup_descubierta 1 265. 265. 10350. 0
## 7 Residuals 32123 821. 0.0256 NA NA
De nuestro modelo confirmamos que todos los coeficientes resultan significativos.
Ahora comparamos nuestro nuevo modelo_log con nuestro modelo_sup_descubierta.
modelos = list(
modelo_log = modelo_log,
modelo_sup_descubierta = modelo_sup_descubierta
)
purrr::map_df(modelos, broom::glance, .id = "model")
## # A tibble: 2 x 13
## model r.squared adj.r.squared sigma statistic p.value df logLik AIC
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 mode… 0.913 0.913 1.60e-1 41902. 0 8 1.33e4 -26605.
## 2 mode… 0.921 0.921 3.98e+4 41717. 0 9 -3.86e5 771799.
## # … with 4 more variables: BIC <dbl>, deviance <dbl>, df.residual <int>,
## # nobs <int>
Observamos que el mismo explica mejor la variabilidad \(Rˆ2\)
Ahora vamos a analizar el cumplimiento de los supuestos del modelo lineal
plot(modelo_log)
- Residuos vs valores predichos: parece existir una cierta estructura en los datos. La linealidad se respeta. Podemos observar un comportamiento heterocedástico aunque parece ser más uniforme que el del modelo_sup_descubierta. Se resaltan 3 observaciones (observaciones 3636,3709,7346). - Normal QQ plot: los extremos no se ajustan a la distribución teórica, tanto inferior izquierdo como superior derecho. Sin embargo parece que se comportan ligeramente mejor que el modelo_sup_descubierta (parecen estár mas suavizados los extremos). - Residual vs leverage:. se pueden apreciar 2 puntos con un leverage bastante alto y el grafico resalta 3 observaciones.
Vamos a generar un nuevo modelo utilizando como variables \(log(rooms)\), \(log(bathrooms)\), \(log(surface_total)\), property_type y barrios y otro similar a nuestro modelo_log reemplazando la variable barrios por la variable l3.
Generamos nuestro modelo modelo_log_surface_total que buscará inferir el resultado del \(log(price)\)
modelo_log_surface_total = lm(
log(price) ~ log(rooms) + log(bathrooms) + log(surface_total) + property_type + barrios,
data = train_ds
)
tidy_mlst <- tidy(modelo_log_surface_total, conf.int = TRUE)
tidy_mlst
## # A tibble: 8 x 7
## term estimate std.error statistic p.value conf.low conf.high
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 8.38 0.0139 603. 0. 8.35 8.40
## 2 log(rooms) 0.00164 0.00275 0.594 5.53e- 1 -0.00376 0.00703
## 3 log(bathrooms) 0.0772 0.00297 26.0 2.52e-147 0.0714 0.0831
## 4 log(surface_total) 0.935 0.00314 298. 0. 0.928 0.941
## 5 property_typeDepart… 0.130 0.00566 23.0 9.67e-116 0.119 0.141
## 6 property_typePH 0.0441 0.00599 7.36 1.90e- 13 0.0323 0.0558
## 7 barriosprecio_bajo -0.732 0.00258 -284. 0. -0.737 -0.727
## 8 barriosprecio_medio -0.362 0.00207 -175. 0. -0.366 -0.358
En nuestro modelo podemos observar como el coeficiente de \(log(rooms)\) parece no ser significativo (su p-valor = 0.5) y su intervalo de confianza incluye el 0. El resto de nuestros coeficientes parecen ser significativos.
Debido a que neustros coeficientes para la variable \(log(rooms)\) parece que no son significativos, vamos a analizar si la interacción de la variable si lo es para nuestro modelo.
tidy(
anova(modelo_log_surface_total)
)
## # A tibble: 6 x 6
## term df sumsq meansq statistic p.value
## <chr> <int> <dbl> <dbl> <dbl> <dbl>
## 1 log(rooms) 1 3724. 3724. 166879. 0
## 2 log(bathrooms) 1 1550. 1550. 69458. 0
## 3 log(surface_total) 1 1332. 1332. 59708. 0
## 4 property_type 2 263. 131. 5889. 0
## 5 barrios 2 1808. 904. 40523. 0
## 6 Residuals 32124 717. 0.0223 NA NA
Y resulta que si lo es.
Analizamos el cumplimiento de los supuestos de nuestro modelo.
plot(modelo_log_surface_total)
Podemos ver como nuestro modelo se comporta de forma muy similar al modelo_log. Sin embargo hay una diferencia importante al analizar Residuos vs Leverage. Aca podemos observar como hay una especie de corte en el leverage entre los valores de 0.0007 y 0.0013 (aproximadamente).
Generamos nuestro modelo modelo_log_surface_total que buscará inferir el resultado del \(log(price)\)
modelo_log_l3 = lm(
log(price) ~ log(rooms) + log(bathrooms) + log(surface_covered) + property_type + sup_descubierta + l3 ,
data = train_ds
)
tidy_mll3 <- tidy(modelo_log_l3, conf.int = TRUE)
tidy_mll3
## # A tibble: 63 x 7
## term estimate std.error statistic p.value conf.low conf.high
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 8.53 0.0265 322. 0. 8.48 8.58
## 2 log(rooms) -0.00924 0.00429 -2.15 3.12e- 2 -0.0177 -0.000834
## 3 log(bathrooms) 0.171 0.00431 39.7 0. 0.163 0.180
## 4 log(surface_covere… 0.783 0.00506 155. 0. 0.773 0.793
## 5 property_typeDepar… 0.198 0.00840 23.6 2.20e-122 0.182 0.215
## 6 property_typePH 0.0476 0.00873 5.46 4.87e- 8 0.0305 0.0648
## 7 sup_descubierta 0.00390 0.0000907 42.9 0. 0.00372 0.00407
## 8 l3Agronomía 0.00514 0.0342 0.151 8.80e- 1 -0.0618 0.0721
## 9 l3Almagro 0.00401 0.0166 0.242 8.09e- 1 -0.0286 0.0366
## 10 l3Balvanera -0.153 0.0176 -8.69 3.91e- 18 -0.188 -0.119
## # … with 53 more rows
En este modelo, podemos observar el mismo problema que surgía al incorporar al analisis la variable l3. Sin embargo, el resto de nuestros coeficientes parecen ser estadisticamente significativos.
Debido a que neustros coeficientes para la variable l3 parece que no son significativos, vamos a analizar si la interacción de la variable si lo es para nuestro modelo.
tidy(
anova(modelo_log_l3)
)
## # A tibble: 7 x 6
## term df sumsq meansq statistic p.value
## <chr> <int> <dbl> <dbl> <dbl> <dbl>
## 1 log(rooms) 1 3724. 3724. 80828. 0
## 2 log(bathrooms) 1 1550. 1550. 33642. 0
## 3 log(surface_covered) 1 1381. 1381. 29982. 0
## 4 property_type 2 172. 86.2 1871. 0
## 5 sup_descubierta 1 81.3 81.3 1765. 0
## 6 l3 56 1008. 18.0 391. 0
## 7 Residuals 32069 1477. 0.0461 NA NA
Y resulta que si lo es.
Analizamos el cumplimiento de los supuestos de nuestro modelo.
plot(modelo_log_l3)
- Residuos vs valores predichos: parece existir una cierta estructura en los datos. La linealidad se respeta. Podemos observar un comportamiento heterocedástico (muy similar al modelo_log). Se resaltan 3 observaciones (observaciones 7346,7942, 28022). - Normal QQ plot: los extremos no se ajustan a la distribución teórica, tanto inferior izquierdo como superior derecho. Sin embargo parece haber cierto suavizado en el extremo superior derecho en comparación a otros modelos. - Residual vs leverage:. se pueden apreciar 3 puntos con un leverage bastante alto y el grafico resalta 3 observaciones (3636, 22246,26579).
Ahora vamos a elegir 2 de los modelos previamente desarrollados para evaluar y seleccionar, junto a nuestros 2 nuevos modelos, cual es el que mejor permite inferir el precio de una propiedad. Vamos a elegir los modelos que mejor expresan la variabilidad \(Rˆ2\). Estos son el modelo_log y modelo_sup_descubierta.
modelos = list(
modelo_log = modelo_log,
modelo_sup_descubierta = modelo_sup_descubierta,
modelo_log_surface_total = modelo_log_surface_total,
modelo_log_l3 = modelo_log_l3
)
purrr::map_df(modelos, broom::glance, .id = "model") %>%
arrange(adj.r.squared)
## # A tibble: 4 x 13
## model r.squared adj.r.squared sigma statistic p.value df logLik AIC
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 mode… 0.843 0.842 2.15e-1 2772. 0 62 3.88e3 -7638.
## 2 mode… 0.913 0.913 1.60e-1 41902. 0 8 1.33e4 -26605.
## 3 mode… 0.921 0.921 3.98e+4 41717. 0 9 -3.86e5 771799.
## 4 mode… 0.924 0.924 1.49e-1 55553. 0 7 1.55e4 -30987.
## # … with 4 more variables: BIC <dbl>, deviance <dbl>, df.residual <int>,
## # nobs <int>
Podemos observar que el modelo que mejor explica la variabilidad es el modelo modelo_log_surface_total. Observamos el valor del \(Rˆ2\) ajustado.
Cargamos el dataset de testing y le aplicamos nuestras transformaciones.
propiedades_test <- read_csv(here("/ds/ar_properties_test.csv"))
## Parsed with column specification:
## cols(
## id = col_character(),
## l3 = col_character(),
## rooms = col_double(),
## bathrooms = col_double(),
## surface_total = col_double(),
## surface_covered = col_double(),
## price = col_double(),
## property_type = col_character()
## )
propiedades_test = propiedades_test %>%
mutate(
precio_m2 = price / surface_total,
barrios = factor(
case_when(
precio_m2 <= quantile(precio_m2)[2] ~ "precio_bajo",
precio_m2 > quantile(precio_m2)[2] & precio_m2 <= quantile(precio_m2)[4] ~ "precio_medio",
precio_m2 > quantile(precio_m2)[4] ~ "precio_alto"
)
),
sup_descubierta = surface_total - surface_covered
)
propiedades_test
## # A tibble: 13,772 x 11
## id l3 rooms bathrooms surface_total surface_covered price
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Afdc… Vele… 3 2 95 69 199900
## 2 ESzy… Nuñez 1 1 44 38 147000
## 3 R7Is… Alma… 1 1 40 37 77000
## 4 MFng… Alma… 1 1 40 37 92943
## 5 j3TP… Pale… 1 1 32 30 125000
## 6 jtR+… Pale… 1 1 45 35 150000
## 7 UOCj… Pale… 1 1 33 32 110000
## 8 3mG1… Parq… 1 1 70 70 225000
## 9 Nmld… Alma… 2 1 48 48 129900
## 10 LmF0… Reti… 2 1 58 58 158000
## # … with 13,762 more rows, and 4 more variables: property_type <chr>,
## # precio_m2 <dbl>, barrios <fct>, sup_descubierta <dbl>
Ahora vamos a evaluar nuestros modelos en training y evaluar sus predicciones en testing. Para esto es necesario hacer un análisis distinto para los modelos que utilizan \(log\) que para el modelo_sup_descubierta.
modelos_log = list(
modelo_log = modelo_log,
modelo_log_l3 = modelo_log_l3,
modelo_log_surface_total = modelo_log_surface_total
)
predicciones_training_log = map(.x = modelos_log, .f = augment)
map_dfr(
.x = predicciones_training_log,
.f = rmse,
truth = exp(`log(price)`),
estimate = exp(.fitted),
.id="modelo"
) %>% arrange(.estimate)
## # A tibble: 3 x 4
## modelo .metric .estimator .estimate
## <chr> <chr> <chr> <dbl>
## 1 modelo_log_surface_total rmse standard 44762.
## 2 modelo_log rmse standard 47444.
## 3 modelo_log_l3 rmse standard 62055.
eval_train_sd = augment(modelo_sup_descubierta)
rmse(
data = eval_train_sd,
truth = price,
estimate = .fitted
)
## # A tibble: 1 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 rmse standard 39757.
Analizando los RMSE de nuestros modelos en training, podemos observar que el modelo con el menor valor es nuestro modelo_sup_descubierta seguido (en orden) por los modelos modelo_log_surface_total, modelo_log, modelo_log_l3.
Ahora vamos a evaluar nuestros modelos en testing.
# Aplicamos la función augment a los 4 modelos con el set de testing
predicciones_testing = map(
.x = modelos_log,
.f = augment,
newdata = propiedades_test
)
# Obtenemos el RMSE para los 4 modelos
map_dfr(
.x = predicciones_testing,
.f = rmse,
truth = price,
estimate = exp(.fitted),
.id="modelo"
) %>%
arrange(.estimate)
## # A tibble: 3 x 4
## modelo .metric .estimator .estimate
## <chr> <chr> <chr> <dbl>
## 1 modelo_log_surface_total rmse standard 43832.
## 2 modelo_log rmse standard 46153.
## 3 modelo_log_l3 rmse standard 60873.
pred_log = augment(
modelo_sup_descubierta,
newdata=propiedades_test
)
pred_log
## # A tibble: 13,772 x 13
## id l3 rooms bathrooms surface_total surface_covered price
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Afdc… Vele… 3 2 95 69 199900
## 2 ESzy… Nuñez 1 1 44 38 147000
## 3 R7Is… Alma… 1 1 40 37 77000
## 4 MFng… Alma… 1 1 40 37 92943
## 5 j3TP… Pale… 1 1 32 30 125000
## 6 jtR+… Pale… 1 1 45 35 150000
## 7 UOCj… Pale… 1 1 33 32 110000
## 8 3mG1… Parq… 1 1 70 70 225000
## 9 Nmld… Alma… 2 1 48 48 129900
## 10 LmF0… Reti… 2 1 58 58 158000
## # … with 13,762 more rows, and 6 more variables: property_type <chr>,
## # precio_m2 <dbl>, barrios <fct>, sup_descubierta <dbl>, .fitted <dbl>,
## # .resid <dbl>
rmse(
data = pred_log,
truth = price,
estimate = .fitted
)
## # A tibble: 1 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 rmse standard 38975.
Volvemos a observar como el modelo con el menor valor de RMSE es nuestro modelo_sup_descubierta.
En lineas generales no podemos afirmar que nuestros modelos cumplen los supuestos del modelo lineal.
Sin embargo, podemos utilizar nuestros modelos para predecir los valores del precio de las propiedades. Ahora bien, para elegir el mejor de los modelos propuestos, vamos a fundamentar nuestra decisión en base a la evaluación de las 2 métricas que estamos observando: - \(Rˆ2\) - RMSE
Los modelos que mejor explican la variabilidad medida por \(Rˆ2\) ajustado son el modelo_log, modelo_sup_descubierta y modelo_log_surface_total. Todos estos modelos se encuentran con valores superiores a 0.91.
Sin embargo, nosotros queremos predecir nuevos datos y para ello es importante medir el RMSE para evaluar el error en la predicción. Al analizar esta metrica, observamos un salto importante entre nuestro modelo_sup_descubierta (38974.87) y sus seguidores modelo_log_surface_total (43832.07) y modelo_log (46152.97).
Por ende seleccionaría como mi mejor modelo el modelo modelo_sup_descubierta.